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Experimental methods based on single particle tracking (SPT) are being increasingly employed 
in the physical and biological sciences, where nanoscale objects are visualized with high temporal 
and spatial resolution. SPT can probe interactions between a particle and its environment but the 
price to be paid is the absence of ensemble averaging and a consequent lack of statistics. Here we 
address the benchmark question of how to accurately extract the diffusion constant of one single 
Brownian trajectory. We analyze a class of estimators based on weighted functionals of the square 
displacement. For a certain choice of the weight function these functionals provide the true ensemble 
averaged diffusion coefficient, with a precision that increases with the trajectory resolution. 

PACS numbers: 05.40. Jc, 31.15.xk, 87.16.dp, 61.43.Er 



Single particle tracking (SPT) generates the time se- 
ries of the position of an individual particle trajectory 
B t in a medium (see, e.g., [H, Q)- Properly interpreted, 
the information so obtained provides an insight into the 
mechanisms and forces that drive or constrain the mo- 
tion of the particle [3| . Nowadays single particle tracking 
is extensively used to characterize the microscopic rhe- 
ological properties of complex media 0] and to probe 
the active motion of biomolecular motors [f| . In biologi- 
cal cells and complex fluids, SPT methods have become 
instrumental in demonstrating deviations from standard 
Brownian motion (BM) [6l4l0| . 

The reliability of the information drawn from SPT 
analysis is not always clear: data is obtained at high 
temporal and spatial resolution but at the expense of 
statistical sample size. Time averaged quantities associ- 
ated with a given trajectory are subject to large fluctu- 
ations across trajectories. For a wide class of anomalous 
diffusion problems, for instance, time-averages of certain 
particle's observables are, by their very nature, random 
variables distinct from their ensemble averages 



Even though standard BM is much better understood 
than anomalous diffusion processes, averaging problems 
persist and complicate the analysis of single trajectories. 
Moreover, in bounded systems, substantial manifesta- 
tions of sample-to-sample fluctuations occur in first pas- 
sage time phenomena [1^|. Standard fitting procedures 
applied to a finite Brownian trajectory unavoidably lead 
to fluctuating estimates Df of the diffusion coefficient, 
due to different thermal histories, particle interactions 
with different defects, or simply due to blur and localiza- 
tion errors, as discussed in [16Hl8| | . In fact, Df might be 
very different from the true ensemble average value D, 
as noticed in SPT measurements of diffusion along DNA 



[l9j] . in the plasma membrane Q or in the cytoplasm of 
mammalian cells [20j . 

The broad dispersion of estimate values extracted from 
common SPT analysis raises an important question: 
Does an optimal methodology able to determine the dif- 
fusion coefficient from just one single trajectory exist? 
Clearly, it is highly desirable to have an estimator of this 
kind even for hypothetical pure cases, such as the uncon- 
strained standard BM with perfectly known location at a 
given time. Such an estimator should possess an ergodic 
property i.e., its most probable value should converge to 
the ensemble average and its variance should vanish as 
the observation time increases. Besides, the knowledge of 
the distribution of a family of estimators could provide a 
way to disentangle the effects of the medium complexity 
or localization errors from variations due to the thermal 
noise driving microscopic diffusion. 

In this paper we study the ergodic properties of 
weighted one-time fits to a BM trajectory. We focus on 
a family of weighted least-squares estimators (u M ) of the 
diffusion coefficient of standard d-dimensional BM, given 
by the following quadratic functionals of a trajectory B t : 



(i) 



where B t=0 = 0, u>(t) is some "trial" weight function of 
the form: 



;(*) = 



1 



(t + ty 



(2) 



/i being a tunable exponent, T - the observation time, to 
- a certain lag time (to >C T) and - the normalization 
constant. The term "least-squares" and the choice of the 
weight function will be made clear below. 
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Here we evaluate the distribution P(u^) for arbitrary \i 
and spatial dimension d. To easily compare the accuracy 
of estimators with different values of fi, we chose such 
that E{u M } = 1, where the symbol E{. . .} denotes the 
ensemble average. Hence, w M = Df/D with D being the 
diffusion constant, 



D = 



E{B?} 

2dt 



(3) 



and Df - its estimate from B t . The best choice of fi 
should produce P(uS) whose maximum u* is the closest 
to the ensemble averaged value 1 and have the smallest 
variance Var(u^). Ultimately, we seek the choice at which 



is ergodic, 



Di 



D, independently of B t as 



e = ta/T^O. 

Before we proceed, two remarks are in order. First, 
note that \i = — 1 corresponds to the standard least 
square estimate (LSE) of the square displacement @, 0, 
[20 . |2~0 ] . The case /i = 1 arises when the unconditional 
probability of observing the whole trajectory B t is max- 
imized (assuming that it is Brownian). It is the so-called 
maximum likelihood estimate (ML E), known to be more 
accurate than the LSE [13, [23 1 . 

We next give a physical interpretation of the estimators 
in Eq. (fT]). Consider a least squares functional: 



F = 



1 



cu(t) dt 



t 



[B 2 t - 2dD f t) 



(4) 



which we generalize by adding a time dependent weight 
function u(t) (the standard choice -LSE- is u)(t) = t). 
The value of Df that minimizes F is: 



D 



9f= I- / dtu(t)BU/ 
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2dD 



dttu(t) (5) 



One recovers Eq. ([T]) by choosing the weight function 
u)(t) — (to + and identifying the denominator with 
A^. Hence minimizes a functional (j4]) and can be 
referred to as a weighted least-squares estimator. 

The moment generating function $((7) of the random 
variable u M , Eq. (TTJ, is defined as 



$(cr) = E{e-<™"} 



(0) 



This function can be calculated using the Feynman-Kac 
formula (see Refs. 22J, |23[ for more details). For /i 7^ 2, 
we find that to leading order in e = t^/T 



$(a) 



r (z/) ( — 

Xi 



Xi 



-, -d/2 



(7) 



for \i < 2, while for [i > 2 it obeys 



$(<r) 
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(8) 



where 1/ - 1/(2- /i), X i = d(2-n)/2, X 2 = d(u-2)/2{»- 
1) and / p (z) is the modified Bessel function [2J|. 

The variance Var(u M ) of -P(w M ) is obtained by differ- 
entiating Eqs. © or (jSJ twice with respect to a. For 
arbitrary /i ^ 2 it is then given explicitly by 



Var(u M ) = j 



(2-/i)/(3-/i), A*<2, 
( M -2)/(2 M -3), M >2. 



(9) 



The consequence of the latter equation is shown in Fig. [T] 
Unexpectedly, the variance can be made arbitrarily small 
at leading order in e by taking fj, gradually closer to 2, 
either from above or from below! The slopes at fi = 2 + 
and fj, = 2~ appear to be the same, so that the accuracy 
of the estimator will be the same for approaching /i = 2 
from above or below. 
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FIG. 1. (color online) The variance of the distribution P(u,j) 
for different values of fi. The curves correspond to Eq. |9j|. 
The symbols correspond to the values obtained from direct 
numerical simulations of 3D random walks for (from light to 
dark) e = 5 x 10~ 5 , 5 x 1CT 6 and 5 x 1CT 7 . 



A word of caution is now in order. Finite-e corrections 
to the result in Eq. (0) are of order of 0(s 2 ~^) for 1 < fi < 
2. Therefore the asymptotic behavior above can be only 
attained when e exp (— 1/(2 — fi)). In other words, 
the variance can be made arbitrarily small by choosing 
fi closer to 2, but only at the expense of increasing the 
experimental resolution (to or T — > 00). 

To confirm our analytical results we simulated random 
walks on a 3d lattice and computed P(u^) using Eq. (fTJ) 
from a large ensemble of trajectories, for different values 
of fi and different resolution e. For \i < 1.5 or fx > 
2.5, the variance computed numerically is well described 
by Eq. © and is independent of e (Fig. [TJ. Near \x — 
2, corrections due to the finite resolution are noticeable, 
but the numerics clearly show that the variance of the 
distribution P(u^) decreases as e — >• 0. 

The large- and small-u asymptotics of P(u^) can be 
deduced directly from Eqs. ([7]) and ©. For /1 < 2 and 
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FIG. 2. (color online) The distribution P(it M ) for different 
fi 7^ 2 in 3D systems. The curves correspond to numerical in- 
version of Eq. J7]) and the symbols to direct numerical simula- 
tions of random walks: from dark to light to fi = — 1 (circles), 
H — 1 (squares), /x = 3/2 (triangles), and /x = 1.95 (stars). 
The numerical values were obtained for e = 10 -5 , except for 
fj, — 1.95 for which we used e = 10 -7 . Recall that fi — —1 
corresponds to LSE and \i = 1 to MLE E2, . In the 

inset we depict the curves corresponding to the inversion of 
Eq. ©: from dark to light ix = 5, 2.5, 2.1 and 2.05. 



iju <C 1, P{Ujj) shows a singular behavior: 



exp 



4 2 



Ml 



(10) 



The asymptotic behavior for /x > 2 can be obtained from 
Eq. (fT0|) by simply replacing \i — > X2- Note that Eq. (fT0|) 
describes a bell-shaped function with a maximal value 
u* — > 1 when /x — > 2 from above or below for arbitrary d. 
Next, for u M 3> 1 and /x < 2, we find 



,d/2-l 



exp 



XiTi- 



(11) 



where 7^1 is the 1st zero of the Bessel function J v {z) 
24( | . Results for /x > 2 follow from Eq. (ITT|) via the re- 
placements xi — > X2 and 7i-i/,i — > 7— i/,i< As /x gradually 
approaches 2, the distribution becomes increasingly nar- 
row: the left tails vanish because of the divergence of the 
factor 1/1 2 — /*| in the exponential, while the right tails 
vanish because |2 — /x|7?.„ 1 and |2 — /i|7i-t/ 1 diverge. 

The distributions P(it^), obtained by inverting 
Eqs. ||7J| and ©, are plotted in Fig. [3] Indeed, the maxi- 
mal value u* — > 1 when /x — » 2 either from above or from 
below. Already for /x = 1.95 (or /x = 2.05) we get the 
most probable value u* 0.94, which outperforms the 
LSE (u* « 0.47) and the MLE (u* « 0.6). For ^ = 1.95 
the variance Var(w^) ps 0.032, which is an order of magni- 
tude less than the variances observed for LSE (= 0.5) and 



the MLE (« 0.33). Similarly to Fig. [TJ finite- resolution 
corrections are negligible for /x < 3/2, and P(u^) is well 
described by Eq. 0. For /x = 1.95 and finite resolution 
e = 10 -7 , we obtain a broader distribution and with a 
smaller u* than the corresponding to Eq. ([7]) for infinite 
resolution. However, note that the most probable value 
of P(mi.95) that we obtain at finite resolution is ps 0.84, 
which outperforms the LSE and MLE for infinite resolu- 
tion. 

We turn next to the case /i = 2 with e = to/T small 
but finite, seeking the variance and the distribution of 
%=2- We consider a slightly more general form for ui(t); 



w(i) 




for t < to, 

for t < t < T, 



(12) 



where £ is a tunable amplitude. For such a choice, the 
moment generating function is given explicitly by 



3(a) = 



-d/2 



(13) 



t>± = (S±i) fch(v^ye 



T ± 



2727^0" 
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where <5 = a/1 + 47<r and 7 
tiating Eq. (fl~3|) . we find 



2/d(£ + ln(l/e)). Differen- 



4 31n(l/ £ )-3(l- £ ) + 2(l- £ K + ^ 2 

Var(W2) = 35 (TWM? 



It is a non-monotonic function of £ with a minimum at 
(2 + e)ln(l/e)-3(l-e) 



pt 



ln(l/e) 



1 



(14) 
at 

(15) 



The corresponding optimized variance is given by: 



Var opt (w 2 ) = 



31n(l/e) -4 + 5e-e 2 



3d ln(l/e) (ln(l/e) + 1 + 2e) - 3(1 - e) ' 

(16) 

From Eq. ([TT3)) we find that in 3d Var op t(w2) ~ 
0.144,0.096,0.082 for e = 10" 3 , 10~ 5 , 10~ 6 , respectively. 
When £ — > 0, Var opt (u2) vanishes as 



Var op t(u 2 ) 



dln(l/e) 



(17) 



Therefore, Var op t(w2) can be made arbitrarily small but 
at expense of a progressively higher resolution. In the 
limit e — > the distribution converges to a delta-function. 
The estimators with zx = 2 are the only, in the family de- 
fined by Eq. ([1]), that possess an ergodic property. This is 
shown in Fig. [3] where we plot P(u 2 ) obtained by numer- 
ically inverting Eq. (fTB"]) for different resolutions. The 
symbols in this figure correspond to numerical simula- 
tions using the weight function of Eq. (p3 
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and 



FIG. 3. (color online) The distribution P{u2) obtained from 
a numerical inversion of Eq. (|13p for 3D systems. The curves 
from the dark to light correspond to e = 10~ 3 , 10~ 4 , 10~ 5 
and 10 -6 . The symbols are results of numerical simulations 
of random walks for e = 1CP 3 (circles) and e = 10~ 4 (squares), 
e = 10~ 5 (triangles) and e = 10~ 6 (diamonds). 



Finally let us consider the case of BM recorded at dis- 
crete time steps — 1, . . . , N, (T = A • N), which is 
an important in its own right problem but also will allow 
us to justify the choice of the weight function in Eq. ([2]). 
We focus on the estimator of a general form 



^ N N 

fi= 2dA£ w i B i-i/£j' w i' (18) 

2=1 2=1 

where u>j now is an arbitrary weight function. The cul- 
minating point of our analysis is to determine, via a vari- 
ational approach, the function uij which yields the lowest 
possible variance Var(u), given from Eq. (|18[) by 

2 N N I N \ 2 

Var(u) = -iJ2^ 5I^femin(fc,j) 2 / ^JUj (19) 

2=1 fc=l \2'=1 / 

where min(fc,j) is the minimum of k and j. Minimizing 

^ N N I N \ 

F = ^ £ w 2 ^w fc min(fc,j) 2 - A ^juij - 1 , 
2=1 fe=i \2=1 / 

(20) 

with respect to each w 3 - (A is a Lagrange multiplier enforc- 
ing the constraint E{u} = 1), we find that the optimal 
weight obeys 



N 

ujj min(fe, j) 2 = A • k , k = 1, . 

2=1 

which can be solved exactly to give 

X = N (J2 



k 



vfe=i 



4fc 2 - 1 



(21) 



(22) 



2A 



LO.j 



N \ 1 

^ k J 4j2 _ 1 



4? 2 - 1 \V iv 

The optimal variance in this case reads 



(23) 



Var(u) = i I £ 



A' 



t (4j 2 - 



(24) 



Therefore, the weight function in Eq. (|23l) minimizes the 
discretized least squares functional in Eq. (I20[) and pro- 
duces an ergodic estimator: the smallest possible vari- 
ance (for the class of estimators defined by Eq. (|T5)) ) 
vanishes as N —¥ oo. Choosing some initial time lag 
and turning to the limit A — > and N — > oo, while 
keeping AA^ = T fixed, the weight function in Eq. (|2"3")l 
converges to the form in Eq. @ with /i = 2, which thus 
justifies our choice of the power-law trial weight func- 
tion for continuous-time Brownian motion. Note that for 
N ^> 1, the leading asymptotic behavior of the variance 
in Eq. ([24| coincides with Eq. ([T7|) . but produces slightly 
higher values of the variance (as the former estimator is 
based on an everywhere discrete process). 

To conclude, we have analyzed the ergodic properties 
and the asymptotic behavior of a family of least-squares 
estimators in Eq. ([1]) . We have demonstrated that the es- 
timators with /i = 2 are the only that possess an ergodic 
property, i.e., they can provide the true ensemble aver- 
aged diffusion coefficient from a single trajectory data 
with any necessary precision, but at expense of a pro- 
gressively higher experimental resolution. Convergence 
to the true ensemble average value appears, however, to 
be rather slow: the variance of such an optimal estima- 
tor vanishes only in proportion to l/ln(T). This means 
that for practical purposes the methods based on two- 
time correlation functions can provide better estimators, 
because the variance of the corresponding estimator de- 
cays faster, as 1/T, even in the presence of localisation 
errors 
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